function [var_r2] = varR2(alphas0,betas0,w,tspan,dt)
% var_r2 = varR2(T,W) simulate decomposed ODEs 100 times and calculate
% var(R2) numerically
%
% Dimension of var_r2,w and tspan
% tspan(1)    tspan(2)  ...  tspan(n-1)  tspan(n)
% w(1)        w(2)      ...  w(n-1)
%             var_r2(1) ...  var_r2(n-2) var_r2(n-1)
%
numSamples = 100;
assert(length(w)==length(tspan)-1,'length(w) should be length(tspan)-1');

if nargin == 2
    dt = 0.01;
end

r2 = zeros(length(w),numSamples);
for j = 1:numSamples
    fprintf(1,'@varR2: #%d sample\n',j);
    alphas = alphas0;
    betas = betas0;
    for i = 1:length(w)
        [r2(i,j),alphas,betas] = solveR2(alphas,betas,w(i),[tspan(i),tspan(i+1)],dt);
    end
end
    
var_r2 = var(r2,ones(1,numSamples),2);
